// QtCore
#include <QtCore/QList>
// Локальные
#include "fem/trianglefem.h"

namespace fem
{
    namespace triangles
    {
        struct Rib
        {
            Rib()
            {
                p1 = p2 = -1;
            }

            Rib(int p1, int p2)
            {
                this->p1 = p1;
                this->p2 = p2;
            }

            int p1;
            int p2;
        };

        bool operator==(const Rib& r1, const Rib& r2)
        {
            return (r1.p1 == r2.p1)&&(r1.p2 == r2.p2);
        }

        //========================================================================================================================
        //========================================================================================================================
        CORE_EXPORT int maximum(const linal::Vector &f)
        {
            Q_ASSERT_X(f.dimension() != 0, "TriangleFEM::maximum", "solution is empty");

            int index(0);
            for(int i = 1; i < f.dimension(); i++)
                if(f[i] > f[index])
                    index = i;
            return index;
        }

        //========================================================================================================================
        //========================================================================================================================
        CORE_EXPORT int minimum(const linal::Vector &f)
        {
            Q_ASSERT_X(f.dimension() != 0, "TriangleFEM::minimum", "solution is empty");

            int index(0);
            for(int i = 1; i < f.dimension(); i++)
                if(f[i] < f[index])
                    index = i;
            return index;
        }

        //========================================================================================================================
        //========================================================================================================================
        CORE_EXPORT geometry::_2d::Point maximum(const TriangleNet &net, const linal::Vector& f)
        {
            Q_ASSERT_X(f.dimension() != 0, "TriangleFEM::maximum", "solution is empty");

            int index(0);
            for(int i = 1; i < f.dimension(); i++)
                if(f[i] > f[index])
                    index = i;
            return net.point(index);
        }

        //========================================================================================================================
        //========================================================================================================================
        CORE_EXPORT geometry::_2d::Point minimum(const TriangleNet &net, const linal::Vector& f)
        {
            Q_ASSERT_X(f.dimension() != 0, "TriangleFEM::minimum", "solution is empty");

            int index(0);
            for(int i = 1; i < f.dimension(); i++)
                if(f[i] < f[index])
                    index = i;
            return net.point(index);
        }

        //========================================================================================================================
        //========================================================================================================================
        CORE_EXPORT void formMatrixProfile(const TriangleNet &net, QVector<int> &ig, QVector<int> &jg, QVector<int> &ig2, QVector<int> &jg2)
        {
            QVector< QList<int> > table(net.numOfPoints());
            QVector< QList<int> > table2(net.numOfPoints());

            //====================================================================================================================
            //заносим в таблицу информацию о том какие узлы связаны с какими поскольку матрица имеет симметричный профиль, то
            //достаточно собрать информацию нужную для построения нижнего треугольника. Т.е. учитываем связь двух узлов если
            //глобальный номер первого узла(a) больше номера второго узла(b). В строке №a хранится информация о всех узлах с
            //которыми связан узел a.
            //====================================================================================================================
            for(int i = 0; i < net.numOfElements(); i++)
                for(int j = 0; j < 3; j++)
                    for(int k = 0; k < 3; k++)
                        if(net.pointIndex(i,j) > net.pointIndex(i,k))
                            if(!table[net.pointIndex(i,j)].contains(net.pointIndex(i,k)))
                                table[net.pointIndex(i,j)].append(net.pointIndex(i,k));

            //====================================================================================================================
            //выделяем память под массив ig
            //====================================================================================================================
            ig.resize(net.numOfPoints() + 1);

            //====================================================================================================================
            //подсчитываем количество элементов в портрете матрицы
            //====================================================================================================================
            int num_of_nonzero_elements(0);
            for(int i = 0; i < net.numOfPoints(); i++)
                num_of_nonzero_elements += table[i].size();

            //====================================================================================================================
            //выделяем память под массив jg
            //====================================================================================================================
            jg.resize(num_of_nonzero_elements);

            //====================================================================================================================
            //формируем ig и jg
            //====================================================================================================================
            ig[0] = 0;
            ig[1] = 0;

            for(int i = 1, j = 0; i < net.numOfPoints(); i++)
            {
                for(int k = 0; k < table[i].size(); k++)
                {
                    //заполняем jg
                    jg[j] = table[i][k];
                    table2[table[i][k]].append(j);
                    j++;
                }

                //заполняем ig
                ig[i + 1] = ig[i] + table[i].size();
            }

            //====================================================================================================================
            //выделяем память под массив ig
            //====================================================================================================================
            ig2.resize(net.numOfPoints() + 1);

            //====================================================================================================================
            //подсчитываем количество элементов в портрете матрицы
            //====================================================================================================================
            int num_of_nonzero_elements2(0);
            for(int i = 0; i < net.numOfPoints(); i++)
                num_of_nonzero_elements2 += table2[i].size();

            //====================================================================================================================
            //выделяем память под массив jg
            //====================================================================================================================
            jg2.resize(num_of_nonzero_elements2);

            //====================================================================================================================
            //формируем ig и jg
            //====================================================================================================================
            ig2[0] = 0;

            for(int i = 0, j = 0; i < net.numOfPoints(); i++)
            {
                for(int k = 0; k < table2[i].size(); k++)
                {
                    //заполняем jg
                    jg2[j] = table2[i][k];
                    j++;
                }

                //заполняем ig
                ig2[i + 1] = ig2[i] + table2[i].size();
            }
        }

        //========================================================================================================================
        //========================================================================================================================
        CORE_EXPORT double functionValue(int triangle, const TriangleNet &net, const linal::Vector &f, const geometry::_2d::Point &point)
        {
            double D[3][3], detD;
            detD = (net.point(triangle, 1).x() - net.point(triangle, 0).x())*(net.point(triangle, 2).y() - net.point(triangle, 0).y()) -
                    (net.point(triangle, 2).x() - net.point(triangle, 0).x())*(net.point(triangle, 1).y() - net.point(triangle, 0).y());

            D[0][0] = (net.point(triangle, 1).x()*net.point(triangle, 2).y() - net.point(triangle, 2).x()*net.point(triangle, 1).y())/detD;
            D[0][1] = (net.point(triangle, 1).y() - net.point(triangle, 2).y())/detD;
            D[0][2] = (net.point(triangle, 2).x() - net.point(triangle, 1).x())/detD;
            D[1][0] = (net.point(triangle, 2).x()*net.point(triangle, 0).y() - net.point(triangle, 0).x()*net.point(triangle, 2).y())/detD;
            D[1][1] = (net.point(triangle, 2).y() - net.point(triangle, 0).y())/detD;
            D[1][2] = (net.point(triangle, 0).x() - net.point(triangle, 2).x())/detD;
            D[2][0] = (net.point(triangle, 0).x()*net.point(triangle, 1).y() - net.point(triangle, 1).x()*net.point(triangle, 0).y())/detD;
            D[2][1] = (net.point(triangle, 0).y() - net.point(triangle, 1).y())/detD;
            D[2][2] = (net.point(triangle, 1).x() - net.point(triangle, 0).x())/detD;

            double result(0);
            for(int i = 0; i < 3; i++)
                result += f[net.pointIndex(triangle, i)]*(D[i][0] + D[i][1]*point.x() + D[i][2]*point.y());

            return result;
        }

        //========================================================================================================================
        //========================================================================================================================
        CORE_EXPORT QList<geometry::_2d::Segment> isoline(const TriangleNet &net, const linal::Vector &f, double value)
        {
            QList<geometry::_2d::Segment> isoline;
            QList<int> segments2ribs;
            QList<int> ribs2segments;
            QList<Rib> ribs;

            //====================================================================================================================
            // Находим все сегменты изолинии.
            //====================================================================================================================
            for(int i = 0; i < net.numOfElements(); i++)
            {
                geometry::_2d::Point *begin(NULL), *end(NULL);
                for(int j = 0; j < 3; j++)
                {
                    //TODO: Вынести проверку в отдельную функцию.
                    if((f[net.pointIndex(i, j%3)] >= value && f[net.pointIndex(i, (j + 1)%3)] < value) ||
                            (f[net.pointIndex(i, j%3)] < value && f[net.pointIndex(i, (j + 1)%3)] >= value))
                    {
                        double a = (value - f[net.pointIndex(i, j%3)])/(f[net.pointIndex(i, (j + 1)%3)] - f[net.pointIndex(i, j%3)]);
                        double x = net.point(i, j%3).x() + a*(net.point(i, (j + 1)%3).x() - net.point(i, j%3).x());
                        double y = net.point(i, j%3).y() + a*(net.point(i, (j + 1)%3).y() - net.point(i, j%3).y());
                        if(begin == NULL)
                            begin = new geometry::_2d::Point(x, y);
                        else if(end == NULL)
                        {
                            end = new geometry::_2d::Point(x,y);
                            break;
                        }
                    }
                    if(f[net.pointIndex(i, j%3)] == value && f[net.pointIndex(i, (j + 1)%3)] == value)
                    {
                        begin = new geometry::_2d::Point(net.point(i, j%3));
                        end = new geometry::_2d::Point(net.point(i, (j + 1)%3));
                        break;
                    }
                }

                if(begin != NULL && end != NULL)
                    isoline.append(geometry::_2d::Segment(*begin, *end));
                delete begin;
                delete end;
            }

            //====================================================================================================================
            // Упорядочиваем сегменты так чтобы конец i-го сегмента совпадал с началом i+1 -го сегмента.
            //====================================================================================================================
            for(int i = 0; i < isoline.size(); i++)
                for(int j = i + 1; j < isoline.size(); j++)
                {
                    if(isoline[i].end() == isoline[j].begin())
                    {
                        isoline.swap(i + 1, j);
                        break;
                    }
                    else if(isoline[i].end() == isoline[j].end())
                    {
                        isoline[j].swap();
                        isoline.swap(i + 1, j);
                        break;
                    }
                }

            return isoline;
        }
    }
}
